--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
      name:  <unnamed>
       log:  C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism_figure.log
  log type:  text
 opened on:   9 Oct 2012, 07:01:16

. #delimit ;
delimiter now ;
. *     ***************************************************************** *;
. *     ***************************************************************** *;
. *       File-Name:      economic_conservatism_figure.do                 *;
. *       Date:           October 1, 2012                                 *;
. *       Author:         MRG                                             *;
. *       Purpose:        Take economic_conservatism.dta and replicate    *;
. *                       the results in Table 2.                         *;
. *           Input File:     economic_conservatism.dta                       *;
. *       Output File:    economic_conservatism_figure.log                *;
. *       Data Output:    none                                            *;
.              *       Previous file:  economic_conservatism.dta                       *;
. *       Machine:        laptop                                                          *;
. *     ****************************************************************  *;
. *     ****************************************************************  *;
. set mem 400m;

Current memory allocation

                    current                                 memory usage
    settable          value     description                 (1M = 1024k)
    --------------------------------------------------------------------
    set maxvar         5000     max. variables allowed           1.947M
    set memory          400M    max. data space                400.000M
    set matsize         400     max. RHS vars in models          1.254M
                                                            -----------
                                                               403.201M

. use "C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta", clear;

. set more off;

. sum;

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
        year |    264370    1994.946     5.85372       1981       2004
       ccode |    264370    383.2626    221.2662          2        920
     country |         0
gov_respon~y |    226573     5.64524     3.05514          1         10
 free_market |     13792    .5100783    .4999165          0          1
-------------+--------------------------------------------------------
 inequality1 |     37839    2.135627    1.171327          1          5
attend_rel~s |    249063    4.335658    2.558971          1          8
        wave |    264370    2.951246    .9936813          1          4
         age |    261392    41.18515    16.31476         15        101
highest_ed~n |    184811    4.433659     2.27252          1          8
-------------+--------------------------------------------------------
social_class |    131988    2.667371    .9835625          1          5
income_scale |    225984    4.682708     2.47316          1         10
income_level |    224694    1.954391    .7923673          1          3
        male |    264175    .4814157    .4996555          0          1
inequality~i |    244467    40.99867    9.505983   17.77497     61.819
-------------+--------------------------------------------------------
      region |    264370    7.419272    2.846142          1         10
attendance~l |    211992    8.404166     6.27391          1         24
attendance~e |    213365      19.842    16.56745          1         80
 sub_saharan |    264370    .0637932    .2443846          0          1
  south_asia |    264370    .0465257    .2106211          0          1
-------------+--------------------------------------------------------
   east_asia |    264370    .0513712    .2207541          0          1
     se_asia |    264370    .0223777    .1479089          0          1
     oceania |    264370    .0169346    .1290267          0          1
 middle_east |    264370    .0865113    .2811182          0          1
latin_amer~a |    264370    .1143662    .3182562          0          1
-------------+--------------------------------------------------------
 east_europe |    264370    .2772743    .4476539          0          1
 west_europe |    264370    .2761319    .4470837          0          1
       wave1 |    264370    .0943261    .2922825          0          1
       wave2 |    264370    .2362863    .4248008          0          1
       wave3 |    264370    .2932027    .4552314          0          1
-------------+--------------------------------------------------------
       wave4 |    264370    .3761849    .4844282          0          1

. desc;

Contains data from C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta
  obs:       264,370                          
 vars:            31                          1 Oct 2012 05:01
 size:    36,483,060 (91.3% of memory free)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
              storage  display     value
variable name   type   format      label      variable label
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
year            int    %8.0g       s020       survey year (s020)
ccode           float  %9.0g                  
country         str48  %48s                   country name
gov_responsib~y byte   %8.0g       e037       government should take less responsibility 1 = yes, 10 = no (reverse of e037)
free_market     byte   %8.0g       e127       free market is right choice 1 = yes, 0 = no (e127)
inequality1     byte   %8.0g       e146       important to reduce inequality 1 = yes, 5 = no (e146)
attend_religi~s byte   %26.0g      relig      how often do you attend religious services (reverse of f028)
wave            byte   %8.0g       s002       survey wave (s002)
age             int    %8.0g       x003       age in years (x003)
highest_educa~n byte   %8.0g       x025       level of education (x025)
social_class    byte   %18.0g      social     self-reported social_class (reverse of x045)
income_scale    byte   %8.0g       x047       income_scale (x047)
income_level    byte   %8.0g       x047r      three levels of income (x047r)
male            float  %9.0g                  female = 0, male = 1
inequality_gini float  %8.0g                  Inequality (gini) by Banones
region          byte   %8.0g                  geographic region (aclp)
attendance_in~l float  %9.0g                  attend_religious_services*income_level
attendance_in~e float  %9.0g                  attend_religious_services*income_scale
sub_saharan     float  %9.0g                  1 = sub-Saharan, 0 otherwise
south_asia      float  %9.0g                  1 = South Asia, 0 otherwise
east_asia       float  %9.0g                  1 = East Asia, 0 otherwise
se_asia         float  %9.0g                  1 = South-East Asia, 0 otherwise
oceania         float  %9.0g                  1 = Oceania, 0 otherwise
middle_east     float  %9.0g                  1 = Middle East, 0 otherwise
latin_america   float  %9.0g                  1 = Latin America, 0 otherwise
east_europe     float  %9.0g                  1 = East Europe, 0 otherwise
west_europe     float  %9.0g                  1 = West Europe, 0 otherwise
wave1           float  %9.0g                  1 = WVS Wave 1, 0 otherwise
wave2           float  %9.0g                  1 = WVS Wave 2, 0 otherwise
wave3           float  %9.0g                  1 = WVS Wave 3, 0 otherwise
wave4           float  %9.0g                  1 = WVS Wave 4, 0 otherwise
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Sorted by:  

. *     ****************************************************************  *;
. *       Create a panel ID variable.                                     *;
. *     ****************************************************************  *;
. egen idn=concat(year ccode);

. encode idn, gen(id);

. *     ****************************************************************  *;
. *     ****************************************************************  *;
. *       Replicate Figure 1a                                             *;
. *     ****************************************************************  *;
. *     ****************************************************************  *;
. gllamm inequality1 attend_religious_services income_level attendance_income_level 
>         male highest_education age, i(id) link(ologit) adapt;

Running adaptive quadrature
Iteration 0:    log likelihood = -41918.993
Iteration 1:    log likelihood = -41801.475
Iteration 2:    log likelihood = -41800.248
Iteration 3:    log likelihood = -41800.179
Iteration 4:    log likelihood = -41798.567
Iteration 5:    log likelihood =  -41797.13
Iteration 6:    log likelihood = -41797.099


Adaptive quadrature has converged, running Newton-Raphson
Iteration 0:   log likelihood = -41797.099  (not concave)
Iteration 1:   log likelihood = -41797.093  
Iteration 2:   log likelihood = -41796.904  
Iteration 3:   log likelihood = -41796.904  
 
number of level 1 units = 31637
number of level 2 units = 31
 
Condition Number = 473.26395
 
gllamm model 
 
log likelihood = -41796.904
 
------------------------------------------------------------------------------
 inequality1 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
inequality1  |
attend_rel~s |   .0666968   .0118476     5.63   0.000     .0434758    .0899177
income_level |   .3487275   .0249256    13.99   0.000     .2998742    .3975808
attendance~l |  -.0315607   .0054166    -5.83   0.000    -.0421771   -.0209444
        male |   .1616671   .0211382     7.65   0.000     .1202369    .2030972
highest_ed~n |   .0795949   .0056416    14.11   0.000     .0685375    .0906522
         age |  -.0096144    .000681   -14.12   0.000    -.0109492   -.0082796
-------------+----------------------------------------------------------------
_cut11       |
       _cons |   .2208871   .1467515     1.51   0.132    -.0667405    .5085147
-------------+----------------------------------------------------------------
_cut12       |
       _cons |   1.485797   .1469893    10.11   0.000     1.197703     1.77389
-------------+----------------------------------------------------------------
_cut13       |
       _cons |   2.832301   .1477372    19.17   0.000     2.542742    3.121861
-------------+----------------------------------------------------------------
_cut14       |
       _cons |   3.904202   .1491484    26.18   0.000     3.611877    4.196528
------------------------------------------------------------------------------
 
 
Variances and covariances of random effects
------------------------------------------------------------------------------

 
***level 2 (id)
 
    var(1): .52292282 (.13397087)
------------------------------------------------------------------------------

 

.         estimates store gllamm2;

. *     ****************************************************************  *;
. *       Figure 1a (top) -- odds ratio graph for religious participation *;
. *       sity across income.                                             *;
. *     ****************************************************************  *;
. estimates restore gllamm2;
(results gllamm2 are active now)

. preserve;

. set seed 10101;

. drawnorm MG_b1-MG_b11, n(10000) means(e(b)) cov(e(V)) clear;
(obs 10000)

. postutil clear;

. postfile mypost odds_hat0 lo0 hi0 using sim, replace;

. noisily display "start";
start

. local a=1;

. while `a'<=3 {;
  2.     {;
  3.         generate x_betahat = MG_b1+MG_b3*(`a');
  4.         gen odds0 = 100*(exp(x_betahat)-1);
  5.         egen oddshat0=mean(odds);
  6.         tempname odds_hat0 lo0 hi0;
  7.         _pctile odds0, p(2.5, 97.5);
  8.     scalar `lo0' = r(r1);
  9.     scalar `hi0' = r(r2);
 10.         scalar `odds_hat0'=oddshat0;
 11.         post mypost (`odds_hat0') (`lo0') (`hi0');
 12.         };
 13.         drop x_betahat odds0 oddshat0;
 14.     local a=`a'+0.01;
 15.     display "." _c;
 16.     };
................................................................................................................................................................................
> .........................
.     postclose mypost;

. restore;

. merge using sim;
(note: you are using old merge syntax; see [D] merge for new syntax)

. summarize _merge;

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
      _merge |    264370    1.001521    .0551261          1          3

. drop _merge;

. gen MV=((_n+99)/100);

. replace  MV=. if _n>201;
(264169 real changes made, 264169 to missing)

. graph twoway hist income_level, percent color(gs12) yaxis(2)
>         ||   line odds_hat0 MV, clwidth(medium) clcolor(black) clpattern(solid)
>         ||   line lo0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||   line hi0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||  ,   
>             xlabel(1 2 3, nogrid labsize(2)) 
>             ylabel(-5 -4 -3 -2 -1 0 1 2 3 4 5, nogrid axis(1) labsize(2))
>             ylabel(0 10 20 30 40, axis(2) nogrid labsize(2))
>             yscale(noline alt)
>             yscale(noline alt axis(2))
>             xscale(noline)
>             legend(off)
>             yline(0, lcolor(black) lwidth(thin))
>             xtitle("", size(2.5)  )
>             ytitle("", axis(2) size(2.5))
>             xsca(titlegap(2))
>             ysca(titlegap(2))
>             scheme(s2mono) graphregion(fcolor(white));

.                         *     ****************************************************************  *;
. *       Figure 1a (bottom) -- odds ratio graph for income across        *;
. *       religious participation.                                        *;
. *     ****************************************************************  *;
.                        drop  odds_hat0 lo0 hi0 MV;

. estimates restore gllamm2;
(results gllamm2 are active now)

. preserve;

. set seed 10101;

. drawnorm MG_b1-MG_b11, n(10000) means(e(b)) cov(e(V)) clear;
(obs 10000)

. postutil clear;

. postfile mypost odds_hat0 lo0 hi0 using sim, replace;

. noisily display "start";
start

. local a=1;

. while `a'<=8 {;
  2.     {;
  3.         generate x_betahat = MG_b2+MG_b3*(`a');
  4.         gen odds0 = 100*(exp(x_betahat)-1);
  5.         egen oddshat0=mean(odds);
  6.         tempname odds_hat0 lo0 hi0;
  7.         _pctile odds0, p(2.5, 97.5);
  8.     scalar `lo0' = r(r1);
  9.     scalar `hi0' = r(r2);
 10.         scalar `odds_hat0'=oddshat0;
 11.         post mypost (`odds_hat0') (`lo0') (`hi0');
 12.         };
 13.         drop x_betahat odds0 oddshat0;
 14.     local a=`a'+0.01;
 15.     display "." _c;
 16.     };
................................................................................................................................................................................
> ..............................................................................................................................................................................
> ..............................................................................................................................................................................
> ..............................................................................................................................................................................
> ...
.     postclose mypost;

. restore;

. merge using sim;
(note: you are using old merge syntax; see [D] merge for new syntax)

. summarize _merge;

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
      _merge |    264370    1.005303    .1028507          1          3

. drop _merge;

. gen MV=((_n+99)/100);

. replace  MV=. if _n>701;
(263669 real changes made, 263669 to missing)

. graph twoway hist attend_religious_services, percent color(gs12) yaxis(2)
>         ||   line odds_hat0 MV, clwidth(medium) clcolor(black) clpattern(solid)
>         ||   line lo0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||   line hi0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||  ,   
>             xlabel(1 2 3 4 5 6 7 8, nogrid labsize(2)) 
>             ylabel(0 5 10 15 20 25 30 35 40 45, nogrid axis(1) labsize(2))
>             ylabel(0 5 10 15 20 25, axis(2) nogrid labsize(2))
>             yscale(noline alt)
>             yscale(noline alt axis(2))
>             xscale(noline)
>             legend(off)
>             yline(0, lcolor(black) lwidth(thin))
>             xtitle("", size(2.5)  )
>             ytitle("", axis(2) size(2.5))
>             xsca(titlegap(2))
>             ysca(titlegap(2))
>             scheme(s2mono) graphregion(fcolor(white));

.             *     ****************************************************************  *;
. *     ****************************************************************  *;
. *       Replicate Figure 1b                                             *;
. *     ****************************************************************  *;
. *     ****************************************************************  *;
. clear;

. clear mata;

. clear matrix;

. use "C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta", clear;

. *     ****************************************************************  *;
. *       Create a panel ID variable.                                     *;
. *     ****************************************************************  *;
. egen idn=concat(year ccode);

. encode idn, gen(id);

. xtreg gov_responsibility attend_religious_services income_level attendance_income_level
>         male highest_education age, i(id) theta;

Random-effects GLS regression                   Number of obs      =    142993
Group variable: id                              Number of groups   =       116

R-sq:  within  = 0.0099                         Obs per group: min =       325
       between = 0.0073                                        avg =    1232.7
       overall = 0.0099                                        max =      4041

Random effects u_i ~ Gaussian                   Wald chi2(6)       =   1428.53
corr(u_i, X)       = 0 (assumed)                Prob > chi2        =    0.0000

------------------- theta --------------------
  min      5%       median        95%      max
0.8251   0.8741     0.9016     0.9335   0.9497

------------------------------------------------------------------------------
gov_respon~y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
attend_rel~s |   .0347845   .0081883     4.25   0.000     .0187358    .0508333
income_level |   .2774759    .019969    13.90   0.000     .2383375    .3166144
attendance~l |  -.0139343    .003809    -3.66   0.000    -.0213998   -.0064688
        male |   .1468475   .0155653     9.43   0.000     .1163401     .177355
highest_ed~n |   .0748199   .0040134    18.64   0.000     .0669538     .082686
         age |  -.0019802   .0005357    -3.70   0.000    -.0030301   -.0009304
       _cons |   4.621252   .0992642    46.56   0.000     4.426697    4.815806
-------------+----------------------------------------------------------------
     sigma_u |  .91166065
     sigma_e |  2.9186703
         rho |  .08889266   (fraction of variance due to u_i)
------------------------------------------------------------------------------

. *     ****************************************************************  *;
. *       Figure 1b (top) -- marginal effect plot for religious           *;
. *       participation across income.                                    *;
. *     ****************************************************************  *;
. generate MV=((_n-1)/100);

. replace  MV=. if _n>302;
(264068 real changes made, 264068 to missing)

. matrix b=e(b);

.  matrix V=e(V);

.  scalar b1=b[1,1];

.  scalar b2=b[1,2];

. scalar b3=b[1,3];

. scalar varb1=V[1,1];

.  scalar varb2=V[2,2];

.  scalar varb3=V[3,3];

. scalar covb1b3=V[1,3];

.  scalar covb2b3=V[2,3];

. scalar list b1 b2 b3 varb1 varb2 varb3 covb1b3 covb2b3;
        b1 =  .03478453
        b2 =  .27747593
        b3 = -.01393426
     varb1 =  .00006705
     varb2 =  .00039876
     varb3 =  .00001451
   covb1b3 = -.00002817
   covb2b3 =  -.0000648

. replace MV=. if _n<101;
(100 real changes made, 100 to missing)

. *     ****************************************************************  *;
. *       Calculate the marginal effect of religious participation.       *;
. *     ****************************************************************  *;
. gen conbx=b1+b3*MV if _n<302;
(264169 missing values generated)

. *     ****************************************************************  *;
. *       Calculate the standard errors for the marginal effect.          *;
. *     ****************************************************************  *;
. gen consx=sqrt(varb1+varb3*(MV^2)+2*covb1b3*MV) if _n<302;
(264169 missing values generated)

.  *     ****************************************************************  *;
. *       Generate upper and lower bounds of the confidence interval.     *;
. *     ****************************************************************  *;
. gen ax=1.96*consx;
(264169 missing values generated)

.  gen upperx=conbx+ax;
(264169 missing values generated)

.  gen lowerx=conbx-ax;
(264169 missing values generated)

. gen yline=0;

. graph twoway hist income_level, percent color(gs12) yaxis(2) 
>         ||   line conbx   MV, clpattern(solid) clwidth(medium) clcolor(blue) clcolor(black) yaxis(1)  
>         ||   line upperx  MV, clpattern(dash) clwidth(thin) clcolor(black) 
>         ||   line lowerx  MV, clpattern(dash) clwidth(thin) clcolor(black) 
>         ||   line yline  MV,  clwidth(thin) clcolor(black) clpattern(solid) 
>         ||   ,  
>              xlabel(1 2 3, nogrid labsize(2))  
>              ylabel(-0.02 -0.01 0 0.01 0.02 0.03, nogrid axis(1) labsize(2)) 
>              ylabel(0 10 20 30 40,  axis(2) nogrid labsize(2)) 
>              yscale(noline alt) 
>              yscale(noline alt axis(2)) 
>              xscale(noline) 
>              legend(off) 
>              yline(0, lcolor(black) lwidth(thin)) 
>              xtitle("" , size(2.5)  ) 
>              ytitle("", axis(2) size(2.5)) 
>              xsca(titlegap(2)) 
>              ysca(titlegap(2)) 
>              scheme(s2mono) graphregion(fcolor(white));

.                                      drop conbx consx ax upperx lowerx MV;

. *     ****************************************************************  *;
. *       Figure 1b (bottom) -- marginal effect plot for income across    *;
. *       religious participation.                                        *;
. *     ****************************************************************  *;
. clear;

. clear mata;

. clear matrix;

. use "C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta", clear;

. *     ****************************************************************  *;
. *       Create a panel ID variable.                                     *;
. *     ****************************************************************  *;
. egen idn=concat(year ccode);

. encode idn, gen(id);

. xtreg gov_responsibility income_level attend_religious_services attendance_income_level 
>         male highest_education age,  i(id) theta;

Random-effects GLS regression                   Number of obs      =    142993
Group variable: id                              Number of groups   =       116

R-sq:  within  = 0.0099                         Obs per group: min =       325
       between = 0.0073                                        avg =    1232.7
       overall = 0.0099                                        max =      4041

Random effects u_i ~ Gaussian                   Wald chi2(6)       =   1428.53
corr(u_i, X)       = 0 (assumed)                Prob > chi2        =    0.0000

------------------- theta --------------------
  min      5%       median        95%      max
0.8251   0.8741     0.9016     0.9335   0.9497

------------------------------------------------------------------------------
gov_respon~y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
income_level |   .2774759    .019969    13.90   0.000     .2383375    .3166144
attend_rel~s |   .0347845   .0081883     4.25   0.000     .0187358    .0508333
attendance~l |  -.0139343    .003809    -3.66   0.000    -.0213998   -.0064688
        male |   .1468475   .0155653     9.43   0.000     .1163401     .177355
highest_ed~n |   .0748199   .0040134    18.64   0.000     .0669538     .082686
         age |  -.0019802   .0005357    -3.70   0.000    -.0030301   -.0009304
       _cons |   4.621252   .0992642    46.56   0.000     4.426697    4.815806
-------------+----------------------------------------------------------------
     sigma_u |  .91166065
     sigma_e |  2.9186703
         rho |  .08889266   (fraction of variance due to u_i)
------------------------------------------------------------------------------

. generate MV=((_n-1)/100);

. replace  MV=. if _n>802;
(263568 real changes made, 263568 to missing)

. matrix b=e(b);

.  matrix V=e(V);

.  scalar b1=b[1,1];

.  scalar b2=b[1,2];

. scalar b3=b[1,3];

. scalar varb1=V[1,1];

.  scalar varb2=V[2,2];

.  scalar varb3=V[3,3];

. scalar covb1b3=V[1,3];

.  scalar covb2b3=V[2,3];

. scalar list b1 b2 b3 varb1 varb2 varb3 covb1b3 covb2b3;
        b1 =  .27747593
        b2 =  .03478453
        b3 = -.01393426
     varb1 =  .00039876
     varb2 =  .00006705
     varb3 =  .00001451
   covb1b3 =  -.0000648
   covb2b3 = -.00002817

. replace MV=. if _n<101;
(100 real changes made, 100 to missing)

. *     ****************************************************************  *;
. *       Calculate the marginal effect of income.                        *;
. *     ****************************************************************  *;
. gen conbx=b1+b3*MV if _n<802;
(263669 missing values generated)

. *     ****************************************************************  *;
. *       Calculate the standard errors for the marginal effect.          *;
. *     ****************************************************************  *;
. gen consx=sqrt(varb1+varb3*(MV^2)+2*covb1b3*MV) if _n<802;
(263669 missing values generated)

.  *     ****************************************************************  *;
. *       Generate upper and lower bounds of the confidence interval.     *;
. *     ****************************************************************  *;
. gen ax=1.96*consx;
(263669 missing values generated)

.  gen upperx=conbx+ax;
(263669 missing values generated)

.  gen lowerx=conbx-ax;
(263669 missing values generated)

. gen yline=0;

. graph twoway hist attend_religious_services, percent color(gs12) yaxis(2) 
>         ||   line conbx   MV, clpattern(solid) clwidth(medium) clcolor(blue) clcolor(black) yaxis(1)  
>         ||   line upperx  MV, clpattern(dash) clwidth(thin) clcolor(black) 
>         ||   line lowerx  MV, clpattern(dash) clwidth(thin) clcolor(black) 
>         ||   line yline  MV,  clwidth(thin) clcolor(black) clpattern(solid) 
>         ||   ,  
>              xlabel(1 2 3 4 5 6 7 8, nogrid labsize(2))  
>              ylabel(0 0.1 0.2 0.3, nogrid axis(1) labsize(2)) 
>              ylabel(0 5 10 15 20 25,  axis(2) nogrid labsize(2)) 
>              yscale(noline alt) 
>              yscale(noline alt axis(2)) 
>              xscale(noline) 
>              legend(off) 
>              yline(0, lcolor(black) lwidth(thin)) 
>              xtitle("" , size(2.5)  ) 
>              ytitle("", axis(2) size(2.5)) 
>              xsca(titlegap(2)) 
>              ysca(titlegap(2)) 
>              scheme(s2mono) graphregion(fcolor(white));

.                          drop conbx consx ax upperx lowerx MV;

. *     ****************************************************************  *;
. *     ****************************************************************  *;
. *       Replicate Figure 1c                                             *;
. *     ****************************************************************  *;
. *     ****************************************************************  *;
. clear;

. clear mata;

. clear matrix;

. use "C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta", clear;

. *     ****************************************************************  *;
. *       Create a panel ID variable.                                     *;
. *     ****************************************************************  *;
. egen idn=concat(year ccode);

. encode idn, gen(id);

. xtlogit free_market attend_religious_services income_level attendance_income_level 
>         male highest_education age,  i(id) quad(30);

Fitting comparison model:

Iteration 0:   log likelihood = -8605.0322  
Iteration 1:   log likelihood = -8353.7745  
Iteration 2:   log likelihood = -8353.4622  
Iteration 3:   log likelihood = -8353.4622  

Fitting full model:

tau =  0.0     log likelihood = -8353.4622
tau =  0.1     log likelihood = -8143.1307
tau =  0.2     log likelihood = -8143.4225

Iteration 0:   log likelihood = -8143.7682  
Iteration 1:   log likelihood = -8139.7322  
Iteration 2:   log likelihood = -8139.6809  
Iteration 3:   log likelihood = -8139.6808  

Random-effects logistic regression              Number of obs      =     12416
Group variable: id                              Number of groups   =        10

Random effects u_i ~ Gaussian                   Obs per group: min =       639
                                                               avg =    1241.6
                                                               max =      1707

                                                Wald chi2(6)       =    414.47
Log likelihood  = -8139.6808                    Prob > chi2        =    0.0000

------------------------------------------------------------------------------
 free_market |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
attend_rel~s |   .1282325    .023191     5.53   0.000     .0827789    .1736861
income_level |   .2501959   .0467444     5.35   0.000     .1585786    .3418132
attendance~l |  -.0413501    .011122    -3.72   0.000    -.0631487   -.0195514
        male |   .2985696   .0381066     7.84   0.000      .223882    .3732572
highest_ed~n |   .0917819   .0100421     9.14   0.000     .0720998    .1114641
         age |  -.0142621   .0012315   -11.58   0.000    -.0166758   -.0118485
       _cons |  -.6148298   .1951551    -3.15   0.002    -.9973268   -.2323328
-------------+----------------------------------------------------------------
    /lnsig2u |  -1.553592   .4582673                      -2.45178   -.6554051
-------------+----------------------------------------------------------------
     sigma_u |    .459877   .1053733                      .2934964    .7205773
         rho |   .0604014   .0260081                      .0255154    .1363135
------------------------------------------------------------------------------
Likelihood-ratio test of rho=0: chibar2(01) =   427.56 Prob >= chibar2 = 0.000

. *     ****************************************************************  *;
. *       Figure 1c (top) -- odds ratio plot for religious participation  *;
. *       across income.                                                  *;
. *     ****************************************************************  *;
. preserve;

. set seed 10101;

. drawnorm MG_b1-MG_b8, n(10000) means(e(b)) cov(e(V)) clear;
(obs 10000)

. postutil clear;

. postfile mypost odds_hat0 lo0 hi0 using sim, replace;

. noisily display "start";
start

. local a=1;

. while `a'<=3 {;
  2.     {;
  3.         generate x_betahat = MG_b1+MG_b3*(`a');
  4.        gen odds0 = 100*(exp(x_betahat)-1);
  5.         egen oddshat0=mean(odds);
  6.         tempname odds_hat0 lo0 hi0;
  7.         _pctile odds0, p(2.5, 97.5);
  8.     scalar `lo0' = r(r1);
  9.     scalar `hi0' = r(r2);
 10.         scalar `odds_hat0'=oddshat0;
 11.         post mypost (`odds_hat0') (`lo0') (`hi0');
 12.         };
 13.         drop x_betahat odds0 oddshat0;
 14.     local a=`a'+0.01;
 15.     display "." _c;
 16.     };
................................................................................................................................................................................
> .........................
.     postclose mypost;

. restore;

. merge using sim;
(note: you are using old merge syntax; see [D] merge for new syntax)

. summarize _merge;

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
      _merge |    264370    1.001521    .0551261          1          3

. drop _merge;

. gen MV=((_n+99)/100);

. replace  MV=. if _n>201;
(264169 real changes made, 264169 to missing)

. graph twoway hist income_level, percent color(gs12) yaxis(2)
>         ||   line odds_hat0 MV, clwidth(medium) clcolor(black) clpattern(solid)
>         ||   line lo0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||   line hi0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||  ,   
>             xlabel(1 2 3, nogrid labsize(2)) 
>             ylabel(-4 -2 0 2 4 6 8 10 12, nogrid axis(1) labsize(2))
>             ylabel(0 10 20 30 40, axis(2) nogrid labsize(2))
>             yscale(noline alt)
>             yscale(noline alt axis(2))
>             xscale(noline)
>             legend(off)
>             yline(0, lcolor(black) lwidth(thin))
>             xtitle("", size(2.5)  )
>             ytitle("", axis(2) size(2.5))
>             xsca(titlegap(2))
>             ysca(titlegap(2))
>             scheme(s2mono) graphregion(fcolor(white));

.                        *     ****************************************************************  *;
. *       Figure 1c (bottom) -- odds ratio plot for income across         *;
. *       religious participation.                                        *;
. *     ****************************************************************  *;
. clear;

. clear mata;

. clear matrix;

. use "C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.dta", clear;

. *     ****************************************************************  *;
. *       Create a panel ID variable.                                     *;
. *     ****************************************************************  *;
. egen idn=concat(year ccode);

. encode idn, gen(id);

. xtlogit free_market income_level attend_religious_services attendance_income_level male highest_education age,  i(id) quad(30);

Fitting comparison model:

Iteration 0:   log likelihood = -8605.0322  
Iteration 1:   log likelihood = -8353.7745  
Iteration 2:   log likelihood = -8353.4622  
Iteration 3:   log likelihood = -8353.4622  

Fitting full model:

tau =  0.0     log likelihood = -8353.4622
tau =  0.1     log likelihood = -8143.1307
tau =  0.2     log likelihood = -8143.4225

Iteration 0:   log likelihood = -8143.7682  
Iteration 1:   log likelihood = -8139.7322  
Iteration 2:   log likelihood = -8139.6809  
Iteration 3:   log likelihood = -8139.6808  

Random-effects logistic regression              Number of obs      =     12416
Group variable: id                              Number of groups   =        10

Random effects u_i ~ Gaussian                   Obs per group: min =       639
                                                               avg =    1241.6
                                                               max =      1707

                                                Wald chi2(6)       =    414.47
Log likelihood  = -8139.6808                    Prob > chi2        =    0.0000

------------------------------------------------------------------------------
 free_market |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
income_level |   .2501959   .0467444     5.35   0.000     .1585786    .3418132
attend_rel~s |   .1282325    .023191     5.53   0.000     .0827789    .1736861
attendance~l |  -.0413501    .011122    -3.72   0.000    -.0631487   -.0195514
        male |   .2985696   .0381066     7.84   0.000      .223882    .3732572
highest_ed~n |   .0917819   .0100421     9.14   0.000     .0720998    .1114641
         age |  -.0142621   .0012315   -11.58   0.000    -.0166758   -.0118485
       _cons |  -.6148298   .1951551    -3.15   0.002    -.9973268   -.2323328
-------------+----------------------------------------------------------------
    /lnsig2u |  -1.553592   .4582673                      -2.45178   -.6554051
-------------+----------------------------------------------------------------
     sigma_u |    .459877   .1053733                      .2934964    .7205773
         rho |   .0604014   .0260081                      .0255154    .1363135
------------------------------------------------------------------------------
Likelihood-ratio test of rho=0: chibar2(01) =   427.56 Prob >= chibar2 = 0.000

. preserve;

. set seed 10101;

. drawnorm MG_b1-MG_b8, n(10000) means(e(b)) cov(e(V)) clear;
(obs 10000)

. postutil clear;

. postfile mypost odds_hat0 lo0 hi0 using sim, replace;

. noisily display "start";
start

. local a=1;

. while `a'<=8 {;
  2.     {;
  3.         generate x_betahat = MG_b1+MG_b3*(`a');
  4.         gen odds0 = 100*(exp(x_betahat)-1);
  5.         egen oddshat0=mean(odds);
  6.         tempname odds_hat0 lo0 hi0;
  7.         _pctile odds0, p(2.5, 97.5);
  8.     scalar `lo0' = r(r1);
  9.     scalar `hi0' = r(r2);
 10.         scalar `odds_hat0'=oddshat0;
 11.         post mypost (`odds_hat0') (`lo0') (`hi0');
 12.         };
 13.         drop x_betahat odds0 oddshat0;
 14.     local a=`a'+0.01;
 15.     display "." _c;
 16.     };
................................................................................................................................................................................
> ..............................................................................................................................................................................
> ..............................................................................................................................................................................
> ..............................................................................................................................................................................
> ...
.     postclose mypost;

. restore;

. merge using sim;
(note: you are using old merge syntax; see [D] merge for new syntax)

. summarize _merge;

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
      _merge |    264370    1.005303    .1028507          1          3

. drop _merge;

. gen MV=((_n+99)/100);

. replace  MV=. if _n>701;
(263669 real changes made, 263669 to missing)

. graph twoway hist attend_religious_services, percent color(gs12) yaxis(2)
>         ||   line odds_hat0 MV, clwidth(medium) clcolor(black) clpattern(solid)
>         ||   line lo0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||   line hi0  MV, clpattern(dash) clwidth(thin) clcolor(black)
>         ||  ,   
>             xlabel(1 2 3 4 5 6 7 8, nogrid labsize(2)) 
>             ylabel(-15 -10 -5 0 5 10 15 20 25 30 35, nogrid axis(1) labsize(2))
>             ylabel(0 5 10 15 20 25, axis(2) nogrid labsize(2))
>             yscale(noline alt)
>             yscale(noline alt axis(2))
>             xscale(noline)
>             legend(off)
>             yline(0, lcolor(black) lwidth(thin))
>             xtitle("", size(2.5)  )
>             ytitle("", axis(2) size(2.5))
>             xsca(titlegap(2))
>             ysca(titlegap(2))
>             scheme(s2mono) graphregion(fcolor(white));

.             *     ****************************************************************  *;
.  *     ****************************************************************  *;
. *       The End                                                         *;
. *     ****************************************************************  *;
.  *     ****************************************************************  *;
. 
end of do-file

. log close
      name:  <unnamed>
       log:  C:\matt\publications\ajps4\replication\data_analysis\economic_conservatism.log
  log type:  text
 closed on:   9 Oct 2012, 08:16:43
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
